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ON Numerical methods for solving the ideal magnetohydrodynamic (MHD) equa- 

tions in more than one space dimension must cither confront the challenge of 
controlling errors in the discrete divergence of the magnetic field, or else be 
faced with nonlinear numerical instabilities. One approach for controlling the 
discrete divergence is through a so-called constrained transport method, which 
t£^ is based on first predicting a magnetic field through a standard finite volume 

solver, and then correcting this field through the appropriate use of a magnetic 

Cj vector potential. In this work we develop a constrained transport method for the 

i i 3D ideal MHD equations that is based on a high-resolution wave propagation 

scheme. Our proposed scheme is the 3D extension of the 2D scheme developed 
by Rossmanith [SI AM J. Sci. Comp. 28, 1766 (2006)], and is based on the 
v ^ high-resolution wave propagation method of Langseth and LeVeque [J. Comp. 

f~^ Phys. 165, 126 (2000)]. In particular, in our extension we take great care to 

\Q maintain the three most important properties of the 2D scheme: (1) all quan- 

OQ tities, including all components of the magnetic field and magnetic potential, 

r — , are treated as cell-centered; (2) we develop a high-resolution wave propagation 

^^ scheme for evolving the magnetic potential; and (3) we develop a wave limiting 

approach that is applied during the vector potential evolution, which controls 
. . unphysical oscillations in the magnetic field. One of the key numerical difficul- 

t> ties that is novel to 3D is that the transport equation that must be solved for 

the magnetic vector potential is only weakly hyperbolic. In presenting our nu- 
merical algorithm we describe how to numerically handle this problem of weak 
C$ hyperbolicity, as well as how to choose an appropriate gauge condition. The 

resulting scheme is applied to several numerical test cases. 

Keywords: magnctohydrodynamics, constrained transport, hyperbolic 
conservation laws, plasma physics, wave propagation algorithms 
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1. Introduction 

The ideal magnetohydrodynamic (MHD) equations are a common model for 
the macroscopic behavior of collisionless plasma [8, 19, 28]. These equations 
model the fluid dynamics of an interacting mixture of positively and negatively 
charged particles, where each species is assumed to behave as a charged ideal gas. 
Under the MHD assumption, this mixture is taken to be quasi-neutral, meaning 
that as one species moves, the other reacts instantaneously This assumption 
allows one to collapse what should be two sets of evolution equation into a single 
set of equations for the total mass, momentum, and energy of the mixture. 
Furthermore, the resulting dynamics are assumed to happen on slow time scales 
compared to the propagation time of light waves, which yields a simplified set 
of Maxwell equations. Finally, the term ideal refers to the fact that we assume 
the ideal Ohm's law: E = B x u, which removes any explicit resistivity and the 
Hall term 2 . 

All of the above described simplifications conspire to turn the original two- 
fluid plasma model into a system that can be viewed as a modified version of 
the compressible Euler equations from gas dynamics. In particular, the ideal 
MHD system can be written as a system of hyperbolic conservation laws, where 
the conserved quantities are mass, momentum, energy, and magnetic field. Fur- 
thermore, this system is equipped, just as the compressible Euler equations are, 
with an entropy inequality that features a convex scalar entropy and a corre- 
sponding entropy flux. Indeed, the scalar entropy, with some help from the fact 
that the magnetic field is divergence- free, can be used to define entropy variables 
in which the MHD system is in symmetric hyperbolic form [5, 18]. 

As has been noted many times in the literature (e.g., Brackbill and Barnes 
[6], Evans and Hawley [14], and Toth [40]), numerical methods for ideal MHD 
must in general satisfy (or at least control) some discrete version of the divergence- 
free condition on the magnetic field: 

V-B = 0. (1) 

Failure to accomplish this gcnerically leads to a nonlinear numerical instability, 
which often leads to negative pressures and/or densities. Starting with the paper 
of Brackbill and Barnes [6] in 1980, several approaches for controlling errors in 
V • B have been proposed. An in-depth review of many of these methods can 
be found in Toth [40]. We very briefly summarize the main points below. 

Projection methods. Projection methods for ideal MHD are based on a predictor- 
corrector approach for the magnetic field. Some standard finite volume 
or finite difference method is used to solve the ideal MHD equations from 
time t n to t n+1 . The approximate magnetic field that is predicted at 
time t n+l is denoted by B*. In general, B* contains both a nontrivial 



2 E is the electric field, B is the magnetic field, and u is the macroscopic velocity of the 
plasma fluid. 



divergence-free subspace and a divergent subspace; one would like to ex- 
tract only the divergence- free part and discard the divergent part. This 
can be accomplished by setting 

B™ +1 := B* - VV>, (2) 

and forcing V •B™ +1 = 0, which results in a Poisson equation for the scalar 
potential ip: 

VV = V • B*. (3) 

Once ip is computed, the divergence-free magnetic field at time t n+1 is 
taken to be (2) (see Toth [40] and Balsara and Kim [3] for further discus- 
sion) . 

This method is attractive for it allows a variety of methods to be used 
in the prediction step, and then only requires one Poisson solve per time- 
step to correct it. The clear disadvantage of this approach is that it 
requires a global elliptic solve on a problem, ideal MHD, that is purely 
hyperbolic. This could be especially computationally inefficient in the case 
of adaptively refined grids. 

The 8-wave formulation. The linearized ideal MHD equations support seven 
propagating plane- wave solutions 3 and a stationary plane- wave solution 4 . 
The stationary plane-wave solution comes directly from the fact that the 
(nonlinear) MHD equations preserve the divergence constraint: 

(V-B) jt = 0, (4) 

where , t denotes the partial derivative with respect to time. A seemingly 
unrelated fact was proved by Godunov [18]: the ideal MHD can only be 
put in symmetric hyperbolic form if one adds to the ideal MHD equations 
a term that is proportional to V • B (in effect adding a term that, at least 
on the continuous level, is zero). As it turns out, this additional "source 
term" not only allows the equations to be put in symmetric hyperbolic 
form, but it also restores Galilean invariance. With the additional term 
the divergence wave is no longer stationary, and instead, it propagates 
with the fluid velocity: 

(V-B) >t + V-(uV-B)=0. (5) 

From the point of view of numerical methods, the difference between a sta- 
tionary and a propagating divergence wave turns out to be very significant. 
Powell [31] and Powell et al. [32] showed that numerical methods applied 
to the MHD equations with the symmetrizing "source term" were much 



3 The seven propagating waves are made up of two fast magnctosonic, two Alfven, two slow 
magnctosonic, and an entropy wave. 

4 This is the so-called divergence wave. 



more stable than the same methods applied to the original MHD equa- 
tions. The modified form of the MHD equations has come be called the 
8-wave formulation, since this form of the equations supports eight prop- 
agating plane wave solutions. Although this approach has been used with 
some success (see Powell et al. [32]), it does have a significant drawback: 
the 8-wave formulation is non- conservative and difficulties with obtaining 
the correct weak solution have been documented in the literature (see for 
example Toth [40]). 

Hyperbolic divergence cleaning methods. This method was introduced by 
Dedncr et al. [12], and is a close cousin to the above described projection 
method. The basic idea is to again solve for the divergence error in the 
magnetic field. Instead of solving an elliptic equation, however, a damped 
hyperbolic equation is prescribed for the divergence error. This does not 
produce an exact divergence-free magnetic field; however, it allows for the 
divergence error to be propagated and damped away from where it origi- 
nated. The main advantage of this method is that it is easy to implement 
and requires no elliptic solve. The main disadvantage is that this approach 
has two tunable parameters: the speed of propagation of the error and the 
rate at which the divergence error is damped. 

Constrained transport methods. The constrained transport (CT) approach 
for ideal MHD was introduced by Evans and Hawley [14]. The method is 
a modification of Yee's method [47] for electromagnetic wave propagation, 
and, at least in its original formulation, introduced staggered magnetic 
and electric fields. Since the introduction of the CT framework, several 
variants and modifications have been introduced, including the work of 
Balsara [2], Balsara and Spicer [4], Dai and Woodward [11], Fey and Tor- 
rilhon [15], Londrillo and Del Zanna [26], Rossmanith [35], Ryu et al. 
[36], and De Sterck [37]. An overview of many of these approaches, as 
well as the introduction of a few more variants, can be found in Toth 
[40]. In particular, Toth [40] showed that a staggered magnetic field is not 
necessary. 

The CT framework, at least several versions of it, can also be viewed as 
a kind of predictor-corrector approach for the magnetic field. Roughly 
speaking, the idea is to again compute all of the conserved quantities with 
a standard finite volume method. From these computed quantities one 
then constructs an approximation to the electric field through the ideal 
Ohm's law. This electric field can then be used to update the magnetic 
vector potential, which in turn, can be used to compute a divergence-free 
magnetic field (see §3 for more details). 

The main advantages of this approach are that (1) there is no elliptic solve 
and (2) there are no free parameters to choose such as in the hyperbolic 
divergence-cleaning technique. 

We focus in this work on developing a constrained transport method for the 
3D ideal MHD equations based on the high-resolution wave propagation scheme 



of LcVcquc [23] and its 3D extension developed by Langseth and LeVeque [21]. 
Our proposed scheme is the 3D extension of the 2D constrained transport scheme 
developed by Rossmanith [35]. 

The wave propagation scheme is an unsplit finite volume method that achieves 
second-order accuracy and optimal stability 5 in higher-dimensions through the 
use of so-called transverse Riemann solvers. These transverse Riemann solvers 
are based partly on the work of Colella [9] on corner transport upwind (CTU) 
methods. It is worth noting that there has been recent work on CTU methods in 
the context of MHD by Gardiner and Stone [16, 17] and Mignone and Tzeferacos 
[27]. Gardiner and Stone [16, 17] developed a constrained transport approach, 
while Mignone and Tzeferacos [27] considered the hyperbolic divergence clean- 
ing method of [12] in the context of the CTU scheme. The method we propose 
in this work is therefore in the same class of methods as that of Gardiner and 
Stone [16, 17] (i.e., unsplit finite volume methods with constrained transport), 
and to a lesser extent Mignone and Tzeferacos [27] (i.e., unsplit methods for 
MHD), although the details of our base scheme and especially the details of 
how the magnetic field is updated differ greatly from these approaches. 

In particular, the approach for updating the magnetic field that we propose 
in this work is based on a generalization of the 2D constrained transport scheme 
of [35] , which is equipped with all of the following features: 

1. All quantities, including all components of the magnetic field and magnetic 
potential, are treated as cell-centered; 

2. The magnetic potential is evolved alongside the eight conserved MHD 
variables via a modified non-conservative high-resolution wave propaga- 
tion scheme; and 

3. Special limiters are applied in the evolution of the magnetic potential, 
which control unphysical oscillations in the magnetic field. 

The scheme developed in this work extends all three of the above features to 
the 3D case. The new challenge that arises in 3D is that the magnetic potential 
is a vector potential (instead of a scalar as in the 2D case). Furthermore, this 
vector potential obeys, at least under gauge condition we advocate in this work, 
a non-conservative weakly hyperbolic evolution equation. We are thus faced with 
two important challenges: 

1. Developing a wave propagation scheme for this non-conservative weakly 
hyperbolic evolution equation; and 

2. Constructing appropriate limiters that again have the effect of limiting the 
magnetic potential in such a way as to produce a non-oscillatory magnetic 
field. 

After briefly reviewing the MHD equations in §2, we describe an overview of 
our proposed method in §3. One important issue that arises from this approach 



5 The wave propagation scheme is optimally stable in the sense that no other fully explicit 
method that uses at most a 3-point stencil in ID, a 9-point stencil in 2D, and a 27-point 
stencil in 3D can achieve larger maximum Courant numbers. 



is the gauge choice; we discuss several possible choices in §4. Under the gauge 
condition that we choose, the 3D transport equation that must be solved for the 
magnetic potential is only weakly hyperbolic. We describe in detail in §5 how 
to numerically handle this difficulty. The resulting scheme is applied to several 
numerical test cases in §6. 

2. Basic equations 

The ideal magnetohydrodynamic (MHD) equations are a classical model 
from plasma physics that describe the macroscopic evolution of a quasi-neutral 
two-fluid plasma system. Under the quasi-neutral assumption, the two-fluid 
equations can be collapsed into a single set of fluid equations for the total mass, 
momentum, and energy of the system. The resulting equations can be written 
in the following form 6 : 
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where p, pu, and £ are the total mass, momentum, and energy densities of the 
plasma system, and B is the magnetic field. The thermal pressure, p, is related 
to the conserved quantities through the ideal gas law: 



P= (7-1) 



£-^||Bj| 2 -^||u|| 2 



(8) 



where 7 = 5/3 is the ideal gas constant. 

The equation for the magnetic field comes from Faraday's law: 

B, t +VxE = 0, (9) 

where the electric field, E, is approximated by Ohm's law for a perfect conductor: 

E = B x u. (10) 

Note that we have used in the above expression a comma followed by a subscript 
as short-hand notation for partial differentiation. This notation is standard in 
many areas of mathematics, most notably in relativity theory (e.g., [30]). We 
will continue to use this notation throughout this paper. 



6 We use boldface letters to denote vectors in physical space (i.e., R. ), and || • || to denote 
the Euclidean norm of vector in the physical space. Vectors in solution space, such as q £ R , 
where q is the vector of conserved variables for the ideal MHD equations: q = (p, pu, £, B) , 
arc not denoted using boldface letters. 



Under the Ohm's law assumption, we can rewrite Faraday's law in the fol- 
lowing divergence form: 



B t + V x (B x u) = B t + V • (uB - Bu) = 0. 



(11) 



Since the electric field is determined entirely from Ohm's law, we do not require 
an evolution equation for it; and thus, the only other piece that we need from 
Maxwell's equations is the divergence-free condition on the magnetic field (7). 
A complete derivation and discussion of MHD system (6)-(7) can be found in 
several standard plasma physics textbooks (e.g., [8, 19, 28]). 

2.1. V • B = is an involution 

We first note that system (6), along with the equation of state (8), pro- 
vides a full set of equations for the time evolution of all eight state variables: 
(p, pu, £ , B). These evolution equations form a hyperbolic system. In particu- 
lar, the eigenvalues of the flux Jacobian in some arbitrary direction n (||n|| = 1) 
can be written as follows: 



A 1,8 = u • n =F cf : fast magnetosonic waves, 

A 2,7 = u • n =F c a : Alfven waves, 

A 3,6 = u • n =F c s : slow magnetosonic waves, 

A = u • n : entropy wave, 



X b 



: divergence wave, 



(12) 
(13) 
(14) 
(15) 
(16) 
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The eigenvalues are well-ordered in the sense that 

A 1 < A 2 < A 3 < A 4 < A 5 < A 6 < A 7 < A 8 . 
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(17) 
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The fast and slow magnetosonic waves are genuinely nonlinear, while the re- 
maining waves are linearly degenerate. Note that the so-called divergence-wave 
has been made to travel at the speed u • n via the 8- wave formulation described 



in §1, thus restoring Galilean invariance [18, 31, 32]. Note that despite the fact 
that we use the eigenvalues and eigenvectors of the 8-wave formulation of the 
MHD equations, we will still solve the MHD equations in conservative form (i.e., 
without the Godunov- Powell "source term"). 

The additional equation (7) is not needed in the time evolution of the con- 
served variables in the following sense: 

If (7) is true at some time t = T ' , then evolution equation (6) guar- 
antees that (7) is true for all time. 

This result follows from taking the divergence of Faraday's law, which yields: 

(V-B) t =0. (22) 

For this reason, (7) should not be regarded as constraint (such as the V • u = 
constraint for the incompressible Navier-Stokes equations), but rather an invo- 
lution [10]. 

2.2. The role of V • B = in numerical discretizations 

Although V • B = is an involution, and therefore has no dynamic impact 
on the evolution of the exact MHD system, the story is more complicated for 
numerical discretizations of ideal MHD. Brackbill and Barnes [6] gave a physical 
explanation as to why V • B = should be satisfied in some appropriate discrete 
sense: 

1/ V ■ B ^ 0, then the magnetic force, 

F = V-|BB-^||B|| 2 l|, (23) 

in the direction of the magnetic field, will not in general vanish: 

F B = ||Bj| 2 (V-B) ^0. (24) 

If this spurious forcing becomes too large, it can lead to numerical 
instabilities (see for example [6, 35, 40]). 

Toth [41] analyzed the magnetic force for central finite difference methods 
and develop an approach with the property that F • B = on the discrete level. 
This analysis makes it clear that simply having V • B = on the discrete level 
is not sufficient to guarantee that F • B = is satisfied. Londrillo and Del 
Zanna [26] came to the same conclusion and proposed a modified approach for 
computing numerical fluxes. However, in [35] it is argued that by satisfying an 
appropriate discrete form of V • B — the error in F • B is sufficiently controlled; 
and therefore, the additional modifications that are proposed in [26] and [41] 
are not found to be necessary. We omit the details here; and instead, refer the 
reader to [26, 35, 41]. 

Another explanation as to why V • B = should not be ignored in numerical 
discretizations of MHD from a slightly different point-of-view, was offered by 



Barth [5]. Barth's explanation is based on the well-known result of Godunov 
[18] that the MHD entropy density, 

U{q) = -plog(pp-i), (25) 

produces a set of entropy variables, U iq , that do not immediately symmetrize 
the ideal MHD equations. Instead, a symmetric hyperbolic form of ideal MHD 
can only be obtained if an additional term that is proportional to the divergence 
of the magnetic field is included in the MHD equations: 

g,t + V-F(g) + X,„V_B =0, where x (q) = ( 7 - 1) ^— . (26) 

ideal MHD additional term 

By looking at how the entropy behaves on the discrete level, Barth [5] was 
able to prove that certain discontinuous Galerkin discretizations of the ideal 
MHD cqutions could be made to be entropy stable (see Tadmor [42]) if the 
discrete magnetic field where made globally divergence-free. The implication of 
this result is that schemes that do not control errors in the divergence of the 
magnetic field run the risk of becoming entropy unstable. 

3. An unstaggered constrained transport framework 

The constrained transport framework advocated in this work is based on 
using the magnetic vector potential, A : 1R + x 1R 3 — *■ R 3 , whose curl yields the 
magnetic field: 

B = V x A. (27) 

A key feature of the proposed method is that the magnetic vector potential is 
updated alongside the conserved variables: q = (p, pu, £, B) T . 

The use of the magnetic potential in computing solutions to the MHD equa- 
tions goes all the way back to Wilson [46] in 1975, who considered relativis- 
tic 2D axisymmetric problems, and Dorfi [13] in 1986, who considered fully 
three-dimensional flow. These approaches did not use modern shock-capturing 
techniques; and therefore, the computed solutions exhibited strong numerical 
diffusion. The first direct 7 use of the magnetic potential in the context of modern 
shock-capturing schemes was due to Londrillo and Del Zanna [25]. Subsequent 
work on using a magnetic potential with shock-capturing schemes includes Dc 
Sterck [37], Londrillo and Del Zanna [26], and Rossmanith [35]. 

The main focus of this work is to extend the 2D constrained transport 
method introduced by Rossmanith [35] to three-dimensions. We begin with 
an outline of our method that lists all of the key steps necessary to completely 
advance the solution from its current state at time t — t n to its new state at 
time t = t n+l =t n + At: 



^Constrained transport methods in the tradition of Evans and Hawlcy [14] don't directly 
use the magnetic potential, but as is discussed in Toth [40], a discrete magnetic potential is 
certainly still in the background in the definitions of the staggered magnetic and electric field 
values. 



Step 0. Start with the current state: (p n , pu n , £ n , B n , A n ). 

Step 1. Update MHD variables via a standard finite volume scheme: 

(p n ,pu n ,£ n ,B n ) => (p n+1 ,pu n+1 ,£*,B*), 

where the energy and magnetic field values, £* and B*, are given a * 
superscript instead of ra + 1 to emphasize that these values will be modified 
by the constrained transport algorithm before the end of the time step. 

Step 2. Define the time-averaged velocity: u n+ 2 = i (u™ + u™ +1 ). 

Step 3. Using the above calculated velocity, u" + 5, solve the magnetic 
potential version of the induction equation 8 : 

A )4 + (V x A) x u™+5 = -W- 

This updates the vector potential: A™ => A™ +1 . 

Step 4. Compute the new magnetic field from the curl of the vector po- 
tential (B = VxA): 



\A 3 ] n+1 - \A 3 ] n+1 \A 2 V +1 - \A 2 V l+1 

r D nn+l _ l^ lij+lk L lij-lk L lijk+1 L iijk-1 /no , 

V B liik ~ ^7. ^TT-, ' [1 *> 



r / 4ii™+ 1 _ r,4ir+ 1 \A 3 ] n+1 -\A 3 V l+1 

r R 2]"+l _ I lijk+1 L^ 1 Jijfe-1 l^- li+ljk I li-ljk (r)Q -> 

[ J ^ ~ 2A~z 2A^ ' [ ' 

[^2]"+! _r / 42i™+ 1 l"/iir +1 — T/l 1 l ri+1 

r R 31»+l _ L li+ljk I li-ljk I lij + lk I lij-lk ,„„•. 

where the bracket notation, [-]" fc , is used to clearly distinguish between 
vector component superscripts and grid and time point superscripts and 
subscripts. 

Step 5. Set the new total energy density value based on one of the fol- 
lowing options 9 : 

Option 1: Conserve total energy: 

£ n+l = £ * 

Option 2: Keep the pressure the same before and after the con- 
strained transport step (this sometimes helps in preventing the pres- 
sure from becoming negative, although it sacrifices energy conserva- 
tion): 

gn +i = £* + I(j| B ™ +1 || 2 -||Bl 2 ). 



8 See §4 for an explanation of this equation. 

9 In this paper we always choose Option 1 in order to exactly conserve energy. 
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The above described algorithm guarantees that at each time step, the fol- 
lowing discrete divergence is identically zero: 

\B 1 Y i+1 - [r 1 !'^ 1 \B 2 \ n+l - [r 2 1" +1 

r -Dii+l L \i+ljk L li-ljk . L Jij+lfc L" Jij-lfc 

[V • B] «* = 2AS + 2A^ 

L -U j fc+1 L Jtjk— 1 n 



2Az 

4. Vector potential equations and gauge conditions 

The key step in the constrained transport framework as outlined in the 
previous section is Step 3, which requires one to solve an evolution equation 
for the magnetic vector potential. One question that immediately arises: what 
should be chosen for the gauge condition? In this section we briefly discuss 
several gauge conditions and their consequences on the evolution of the magnetic 
vector potential. 

The starting point for this discussion is the induction equation: 

B, 4 +V x (B x u) =0, (32) 

where, for the purposes of the algorithm outlined in §3, we take the velocity, u, 
as a given function. We set B = V x A and rewrite (32) as 

V x {A, t + (V x A) x u} = 0, (33) 

=*> A t |(VxA)xu= -VV>, (34) 

where ip is an arbitrary scalar function. Different choices of tp represent different 
gauge condition choices. 

Before we explore various gauge conditions, however, it is worth pointing 
out that the situation in the pure two-dimensional case (e.g., in the xy-plane) 
is much simpler. The only component of the magnetic vector potential that 
influences the evolution in this case is A 3 (i.e., the component of the potential 
that is perpendicular to the evolution plane); and furthermore, all gauge choices 
lead to the same equation: 

A 3 t + u 1 A% + u 2 A 3 y = 0, (35) 

where A 3 is uniquely defined up to an additive constant. 

4-1. Coulomb gauge 

An obvious choice for the gauge is to take the vector magnetic potential to 
be solenoidal: 

V-A = 0, (36) 
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resulting in the Coulomb gauge. We are now able to add — u V • A to the left- 
hand side of (34) and obtain an equation for the potential that is in symmetric 
hyperbolic form: 
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(37) 



The main difficulty with this approach, however, is that at each time step one 
must solve a Poisson equation to determine the Lagrange multiplier ip: 



-V 2 V> = V • [(V x A) x u] . 



(38) 



Having to solve an elliptic equation in each time step makes this approach have 
the same efficiency problems as the projection method. 

4-2. Lorentz-like gauge 

In the ideal MHD setting, since the speed of light is taken to be infinite, 
the Lorentz and Coulomb gauges are equivalent. However, one possibility is to 
introduce a fictitious wave speed, £, that is larger than all other wave speeds in 
the MHD system. We can then take 



^*=£ 2 V-A, 
which results in the following evolution equation for (A,-0): 
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(40) 



The flux Jacobian of this system in some direction n (where ||n|| = 1) can be 
written as 
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x=i- 


4. 


£, u • n, u 


n). 



rv' 




(41) 



(42) 
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If £ > |u • n|, the right-eigenvectors of N(n) are complete, and thus, the system 
is hyperbolic. 

Although this seems like a potentially useful gauge choice, we found in prac- 
tice that numerical solutions to system (40) did not produce accurate magnetic 
fields. In particular, we observed errors in the location of strong shocks, which 
is presumably due to the fact that on the discrete level 

V x W + 0, (43) 

thus resulting in errors in B. 

4-. 3. Helicity-inspired gauge 

Another gauge possibility is to directly set the scalar function, ip, to some- 
thing useful. One such choice is 

ip = u • A, (44) 

which is inspired by the magnetic helicity: B • A. This yields the system: 

Aj + a'A^I^Ajl^A^MA, (45) 

where 



u 2 A 


y+U 


3 A 


~ u ]x 


u 2 x 


u%~ 


u )y 


u % 


u % 


k 


< 


u %. 



M := - u) y u 2 y u A y . (46) 

u z u z u z 

One obvious approach for solving this equation is via operator splitting, whereby 
equation (45) is split into three decoupled advection equations: 

A it + v}A >x + u 2 A tV + u 3 A. z = 0, (47) 

and a 'linear' ordinary differential equation 10 : 

A )t = MA. (48) 

The main difficulty with this approach is that the matrix M in equation (48) 
could (and often does) have eigenvalues that have a positive real part; thereby, 
causing this system to be inherently unstable. For this reason, numerical tests 
using this gauge condition were generally not successful. 

4-4- Weyl gauge 

The choice that we finally settled on was the Weyl gauge: 

^ = 0. (49) 



10 This equation is linear in the sense that the velocity, u, is taken to be frozen in time at 
i = t „+l/2_ 
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In this approach, the resulting evolution equation is simply (34) with a zero 
right-hand side. This gauge is the most commonly used one in the description 
of constrained transport methods (see for example [25, 40]). 

As we will describe in detail in the next section, §5, the resulting system is 
only weakly hyperbolic. This is due to the fact that there are certain directions 
in which the matrix of right-eigenvectors of the flux Jacobian does not have full 
rank. This degeneracy causes some numerical difficulties, which we were able to 
overcome through the creation of a modified wave propagation scheme 11 . The 
details of this approach are described in the next section. 

5. Numerical methods 

The numerical methods developed in this paper are based on the high- 
resolution wave propagation scheme of LeVeque [23] and its 3D extension de- 
veloped by Langseth and LeVeque [21]. In §5.1 we briefly review the 3D wave 
propagation approach. In §5.2 we show in detail how this approach can be mod- 
ified to solve the non-conservative and only weakly hyperbolic magnetic vector 
potential equation (62)-(63). In §5.3 we develop a limiting strategy that is ap- 
plied during the magnetic potential update, but designed to control unphysical 
oscillations in the magnetic field. Finally, in §5.4 we briefly mention how our 
approach simplifies in the special case of 2.5-dimensional problems. 

5.1. The wave propagation scheme of Langseth and LeVeque [21] 

In Step 1 of the constrained transport method described in §3 we apply a 
numerical method for the three-dimensional MHD equations. Here we use a ver- 
sion of the three-dimensional wave propagation algorithm of Langseth and LeV- 
eque [21], see also [24], which is based on a decomposition of flux differences at 
grid cell interfaces as outlined in [1] . This is a multidimensional high- resolution 
finite volume method that is second order accurate for smooth solutions and 
that leads to an accurate capturing of shock waves. 

To outline the main steps of this algorithm, we consider a three-dimensional 
hyperbolic system of the general form 

q,t + f{q).x + g(q), y + %),, = 0, (50) 

with q : R + x R 3 — *• R m and f,g,h : R m —$■ W n . In quasilincar form the system 
reads 

q tt + A(q) q tX + B(q) q tV + C{q) q. z = 0, (51) 

where A, B, and C are the flux-Jacobians in each of the three coordinate direc- 
tions: 

A(q) := f(q), q , B(q) := g(q), q , and C(q) := h(q), q . (52) 



11 Sec also Fey and Torrilhon [15] for a discussion of numerical discretizations of the weakly 
hyperbolic induction equation. 
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We construct a 3D Cartesian mesh with grid spacing Ax, Ay, and Az in each 
of the coordinate directions. We represent the solution at each discrete time 
level t = t n as piecewise constant, such that in each grid cell (i,j,k) the solu- 
tion is given by Q" fc - Q"- fc denotes the cell average of the conserved quantity 
q(t n , x, y, z) in the grid cell centered at (cc;, yj, z^): 

1 /■*«+¥ rvj+¥ r^+¥ 



Q7jk ^ a a A / / / q(t n ,Z,r,,0 <%dr,d{. (53) 

J AxAyAz J x ._^ J Vi -^ Jz k -^ 

The numerical method can be written in the form 



Qm = Qlk ^ (A + AQ t _ yk + A-AQ lH jk 
-^(B + A Qij _ hk + B-AQ ij+i 

-^ z ( C+AQ ^ k - l z +C ~ AQ ^ k+ i) (M) 

~ A^ \ F i+hok~ F i-hJk) - -^ \ G n+lk- G ij-ik 
At 



The first three lines in (54) describe a first order accurate update and the last 
two lines represent higher order correction terms. The wave propagation method 
is a so-called truly multidimensional scheme in the sense that no dimensional 
splitting is used to approximate the mixed derivative terms that are required in 
a second order accurate update. 

At each grid cell interface in the x-direction, we decompose the flux differ- 
ences 

AF i _i jk = f(Q? jk )-f(Q?_ ljk ), (55) 

into waves 



Z p 

i-\jk 



e>_ hk -A Fl _ hok 



r 



v 



(56) 



which are moving with speeds s p ± , p = 1, . . . , M w . For this decomposition 

we use the left and right eigenvectors proposed by Roe and Balsara [33] (see also 
Powell [32] for a discussion of these eigenvectors) for a linearized flux Jacobian 
matrix of the MHD equations. The fluctuations used in the first order update 
have the form: 

A + AQ l _ hk = 
A-AQ i _, ik = 



y z p !,+- y z p t . h , 

p:s F , >0 " p:s" . =0 


(57) 


v z p + - y z p 

p:s F , <0 " p:s p , =0 


(58) 
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Analogously, we compute fluctuations in the y— and z— directions. An advan- 
tage of the decomposition of the flux differences compared to the standard wave 
propagation method from [21, 23], which is based on an eigenvector decomposi- 
tion of the jumps of the conserved quantities at each grid cell interface, is that 
the local linearization of the flux Jacobian matrix does not require Roe averages, 
which are quite difficult to compute for the MHD equations [7]. Here we use 
simple arithmetic averaging instead. See Wesenberg [45] for a discussion of how 
the use of arithmetic averages is equivalent to an MHD Riemann solver based 
on van Leer flux- vector splitting [44] . 

The waves and speeds are also used to compute high-resolution correction 
terms. For the waves in the x-dircction, the correction fluxes have the form 



\jk 



= -^sign 



i—hjk 



1 - 



At 

Ax 



i—hjk 



i-h]k^ 






(59) 



where 



i-njk 



Z v ■ Z v 

% 2 3 k i 2 3 k 

~Z? . 2, p ' 



where / = 




if 
if 



i-\jk 

S P 
i-\jk 



>o, 
<o, 



(60) 



is a measure of the location smoothness in the p-th wave family and 
total variation diminishing (TVD) limiter function [20, 43, 39]. 

The wave propagation method is a 1-step, unsplit, second-order accurate 
method in both space and time (for smooth solutions); and therefore, is based 
on a Taylor series expansions in time: 



-(•%/, 



q(t + At,-K.) =q 
q - At f yX + g iV + h tZ 
V{Ah, z ) x + {Bf iX ) v + 



1 

+ 2- 

(Bh, x ) 



At 2 q ttt 



-At 



{W,x) lX 
+ (C/,x) z ■ 



Transverse terms 



f O (At 3 ) 
(Bg, y ) y + (Ch z ) z 

(Cg, y )A + o(At 3 ), 

(61) 



where time derivatives have been replaced by spatial derivations through the 
conservation law. The wave propagation method as described so far, approxi- 
mates each of the terms in the above Taylor series to second order accuracy, ex- 
cept those terms marked as transverse terms. The transverse terms are approx- 
imated via additional Riemann solvers known as transverse Riemann solvers. 
Additionally, Langseth and LeVeque [21] found that in order to achieve opti- 
mal stability bounds, additional so-called double transverse Riemann solvers 
are required to approximate certain mixed third-derivative terms. We omit a 
full discussion of the transverse and double transverse terms in this work, and 
instead, refer the reader to the paper of Langseth and LeVeque [21]. 
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5.2. The evolution of the magnetic potential 

We now describe a numerical method for the evolution equation of the mag- 
netic potential using the Weyl gauge, i.e., for the discretization of (34) with zero 
right hand side. The equation can be written in the form 



A, t + JVi(u) A tX + N 2 (u) A „ + N 3 (u) A z = 0, 



(62) 



with 


























"o 


-u 2 


-u 3 ' 




' u 2 





" 




~u 3 





0" 


m = 





u 1 





,N 2 = 


-u 1 





— u 3 


,N 3 = 





u 3 













u 1 _ 










u 2 _ 




-u 1 


-u 2 






(63) 



5.2.1. Weak hyperbolicity 

We recall the following definitions regarding hyperbolicity. 

Definition. (Strict hyperbolicity) The quasilinear system, 
q yt + A{q) q iX + B(q) q yV + C(q) q. z = 0, 

is strictly hyperbolic if the matrix 

M(n, q) := n 1 A{q) + n 2 B(q) + n 3 C(q) 



(64) 



(65) 



is diagonalizable with distinct real eigenvalues for all ||n|| = 1. In other words, 
the p th eigenvalue of M(n,q) is real and has geometric and algebraic multiplic- 
ities of exactly one. 

Definition. (Non-strict hyperbolicity) The quasilinear system (64) is non- 
strictly hyperbolic if the matrix (65) is diagonalizable with real but not necessarily 
distinct eigenvalues for all ||n| = 1. In other words, the p th eigenvalue of 
M(n, q) is real and has the same geometric and algebraic multiplicity, but that 
multiplicity may be greater than one. 



Definition. (Weak hyperbolicity) The quasilinear system (64) is weakly 
hyperbolic if the matrix (65) has real eigenvalues for all ||n|| = 1, but is not 
necessarily diagonalizable. In other words, the p eigenvalue ofM(n,q) is real 
but may have an algebraic multiplicity larger than its geometric multiplicity. 



The flux Jacobian of system (62) in an arbitrary direction n £ S 2 is 



n x N x +n 2 N 2 +n 3 N 3 



n 2 u 2 + n 3 u 3 



i i 
-n u 



n u + n 3 u 



-n 3 u 2 



The eigenvalues of this system are 

A = {0,n • u,n • u}; 



— n 2 u 3 
n u + n 2 u 2 



(66) 



(67) 
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and therefore, we always have real eigenvalues. The eigenvectors can be written 
in the following form: 



R = 


r (i) 


r (2) 


r (3) 


= 


n 1 

o 

rr 
n 3 


n 2 u 3 — n 3 u 2 u 1 (u • n) 
n 3 U l — n l u 3 u 2 (u • n) 
n^u 2 — r^u 1 u 3 (u • n) 


ssuming that ||u|| ^ and 


|n|| = 1, the determinant of 


ritten 


as 






det 


(R) 


= — ||u|| 3 cos(a) sin 2 (a), 



n 1 


|u 


2" 


n 2 


|u 


2 


n 3 


|u 


2 



(69) 

where a is the angle between the vectors n and u. The difficulty is that for any 
non-zero velocity vector, u, one can always find four directions, a = 0, 7r/2, 7r, 
and 37r/2, such that det(i?) = 0. In other words, for every ||u|| ^ there exists 
four degenerate directions in which the eigenvectors are incomplete. Therefore, 
system (62)-(63) is only weakly hyperbolic. 

5.2.2. An example: the difficulty with weakly hyperbolic systems 

Non-strict hypcrbolicity is common in many standard equations (e.g., Eulcr 
equations of gas dynamics); and, although it causes some difficulties in proving 
long-time existence results, it generally does not cause problems for numerical 
discretization of such equations. Weak hypcrbolicity, however, is a different 
story. We illustrate this point with the following simple example. Let 



,* 







= 0, 



(70) 



where e e R is a constant. The eigen-decomposition of the flux Jacobian can 
be written as follows: 



-e 1 
e 



= RAR- 1 = 



"1 


1" 




— e 


If 


1 


~2e 


-1" 





2c 







£ 


' 2e 





f 



(71) 



Since the eigenvalues are always real, this system is hyperbolic. For all e =/= 0, 
the system is strongly hyperbolic, and for e = the system is only weakly 
hyperbolic. Since we have the exact eigen-decomposition, we can write down 
the exact solution for the Cauchy problem for all e: 



u (x + et) - ^{v (x + et) - v (x — st)} 
vq(x — et) 



In the weakly hyperbolic limit we obtain: 



lim 

e->0 



u (x) -tv' (x) 
v (x) 



(72) 



(73) 



The change in dynamics between the strongly and weakly hyperbolic regimes 
is quite dramatic: in the strongly hyperbolic case the total energy is conserved, 
while in the weakly hyperbolic case the amplitude of the solution grows linearly 
in time. 
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5.2.3. An operator split approach 

In the case of the magnetic vector potential equation (62)-(63), the weak 
hyperbolicity is only an artifact of how we are evolving the magnetic potential. 
One should remember that the underlying equations (i.e., the ideal MHD sys- 
tem) are in fact hyperbolic. One way to understand all of this is to remember 
that the true velocity field, u, depends on the magnetic potential in a nonlinear 
way. That is, even though there might be instantaneous growth in the vector 
potential due to the weak hyperbolicity in the potential equation, this growth 
will immediately change the velocity field in such a way that the overall system 
remains hyperbolic. 

Numerically, however, we must construct a modified finite volume method 
that can handle this short-lived weak hyperbolicity. Through some experimen- 
tation with various ideas, we found that simple operator splitting techniques 
lead to robust numerical methods. In particular, we found that there are two 
obvious ways to split system (34) into sub-problems. 

The first possibility is based on a decomposition of the form 



Sub-problem 1: 



Sub-problem 2: 



Sub-problem 3: 



]t + ^A) y 


+ v?A\ 


= 0, 


-,t " u Ay 


= o, 




-t - «MJ, 


= 0, 




,t — u A» 


= o, 




2 t+u l A 2 x 


+ u 3 A% 


= 0, 


3 t - u 2 A 2 z 


= 0, 




,t ~ u A x 


= o, 




,t ~ u Ay 


= 0, 




3 t + ^A% 


+ u 2 A% 


= 0. 



(74) 



(75) 



(76) 



First, consider Sub-problem 1. The first equation can be approximated using the 
two-dimensional method described in [35, Section 5.3], which is a modification 
of LeVeque's wave propagation method that leads to a second order accurate ap- 
proximation for smooth solutions, both in the potential and the magnetic field, 
and a non-oscillatory solution simultaneously for both the magnetic potential 
and magnetic field. Once the quantity A 1 is updated, central finite difference 
approximations can be used to update A 2 and A 3 . Sub- problems 2 and 3 can be 
approached analogously and Strang operator splitting can be used to approxi- 
mate (62). We have tested this splitting and obtained very satisfying results. 
A second approach to split (34) is based on dimensional splitting, i.e., we 
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consecutively solve the problems 
Sub-problem 1: 



A\ - u 2 A% 


u 3 A% 


= 


A 2 t +u'A% 


= 0, 




A\ + u x A\ 


= 0, 





(77) 



Sub-problem 2: A\ + u 2 A\ = 0, 



Sub-problem 3: 



A l~^ A )y 


« 8 ^ 


= 0, 


4 + " 2 4, 


= 0, 




A] t + u 3 A] z 


= 0, 




A 2 t + u 3 A 2 z 


= o, 




4 - u x A\ 


- u 2 A% 


= 0. 



(78) 



(79) 



We will now present a numerical method that is based on this second splitting. 
Our method should be easy to adapt by other users; and, in particular, also 
by those who do not base their numerical method on the wave propagation 
algorithm. 

We denote the numerical solution operator for Sub-problem 1, 2 and 3 by 
L\, L2 and L3, respectively. Then a Strang- type operator splitting method can 
be written in the form 

A „+l = L M,2 L ft/2 L At L At/2 L At/2 An ^ {gQ) 

One time step from t n to t n+l consists in consecutively solving Sub-problem 1 
over a half time step, Sub-problem 2 over a half time step, Sub-problem 3 over 
a full time step, Sub-problem 2 over another half time step and finally solving 
Sub-problem 1 over another half time step. 

5.2.4- Discretization of Sub-problem 1 

We now present the details for the approximation of the first Sub-problem 
(77) for a time step from r™ to r n+ . The evolution of A 2 and A 3 is described 
by two decoupled scalar transport equations, which have the form 

A* t + u 1 {T n+ K*)A% = 0, (81) 

4 + U 1 (t"+5,x)A 3 x -0. (82) 

We assume that cell average values of A 2 and A 3 are given at time t™ and we 
wish to approximate A 2 and A 3 at time r™ +1 . In our application, cell centered 
values of the advection speed u 1 are given at the intermediate time level (see 
Step 2 of the algorithm from §3). 
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The update for all i,j, k and I <G {2, 3} has the form 



Wm = W 



ijk 



At 

Ax 
At 
Ax 



A~ A A 



j+j j k 



A + AA e 



hik 



f: 



e+ 



i - F: 

-\- 2 J K I o 3 "* 



with 



A~AA{ 

A+AA e 






; j k 



:= \A e ] n - \A e ] n 

I \ ijk L J i — 1 j k ' 

:=min(ul_ ljk ,0)W^i jk , 
:=max(uJ ifc ,0)Vl>f_i J . fc) 



(83) 



(84) 
(85) 
(86) 



and 



F. 



i+\jk 



F 1 



i j k 



l ijk\ 



l ijk\ 



1- 



A7 



i+ j j k 



Aa 



1- 



A7 



*ijfcl 



■jifc 



Ax \u} jk \ 



ijk 



Hjk 



w. 



»+ 5 j k 



m 



Here w; 



5(«i-lrt + 



0(^- fe ), (87) 

i jfe ^(0- (88) 

u i7fe) * s a ^ ace centered 



iyj. is a cell centered value and w*_i . fc — 5 \ u i-ijk~ ru 'ijkJ 
value. The key innovation in this construction as developed in [35] is the use of 
a cell-based limiter, rather than the standard edge-based limitcr as represented 
in equations (59). Furthermore, in this construction the smoothness indicator, 
#|- fe , is based on wave-differences rather than waves as in (60): 



e t = AW^ 



J ijk 



AW. 



(89) 



ijk 



where AWf, t = W £ , 1 ., — W. 1 ., , and the index / is again chosen from the 

upwind direction. The reason for the use of wave-differences instead of waves is 
that ultimately, the physical quantity of importance is the magnetic field and 
not the magnetic potential. As was shown in [35], standard wave limiters do 
not control oscillations in derivative quantities (i.e., the magnetic field), while 
the wave-difference limiter approach does. 

Once we have updated A 2 and A 3 , we update A 1 by discretizing the equation 



A\- 



u 2 A\ - u 3 A 3 = 0. 



For this we use a finite difference method of the form 



l ijk 



D X [A 



21 "+ 1 
ijk 



+h;; 5 (^[ a " 



i]k 



*>.W£ 



(90) 



(91) 
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where D x is the standard centered finite difference operator in the x-direction, 
i.e., 



\A2] n _ \A2V 

,2]™ L ii+ljk I li-ljk 

^' fc '" 2Ax 



^ [A2] n ;= L Ji+lj* L Jj-lj^ (g2) 



We will refer to the numerical solution of (81) and (82) as hyperbolic solves and 
the numerical solution of equation (90) as the weakly hyperbolic solve. 

5.3. Additional limiting in the weakly hyperbolic solve 

For the spatially two-dimensional case, the limiting of wave-differences elim- 
inates unphysical oscillations in the magnetic field. In the three-dimensional 
case, the equation for the magnetic potential is more challenging and only weakly 
hyperbolic. Each sub-problem of the numerical method is a combination of one- 
dimensional hyperbolic solves (i.e., equations (81) and (82)) and a part that we 
refer to as the weakly hyperbolic solve (i.e., equation (90)). The wave-difference 
limiting strategy as outlined in §5.2.4 is generally insufficient in controlling un- 
physical oscillations in the magnetic field variables. In particular, the step that 
produces unphysical oscillations is the weakly hyperbolic solve as described in 
equation (91). Since this part of the update is quite different from standard 
upwind numerical methods for hyperbolic equations, there is no obvious place 
to introduce wave-difference limitcrs. Instead, we describe in this section an 
approach based on adding artificial diffusion to this part of the update in order 
to remove these unphysical oscillations. 

We introduce a diffusive limiter that is inspired by the artificial viscosity 
method that is often used in other numerical schemes such as the discontinuous 
Galerkin approach (e.g., see [29]). Instead of (90), we discretize the problem 
with an added dissipative term 

A] t - u 2 A% - u 3 A% - e l A] xx , (93) 

where the parameter e 1 depends on the solution structure and controls the 
amount of artificial diffusion. We choose 



£ 



1 _ o,,l„l 



Ax 



2 



2-V — , (94) 



where v 1 is a positive constant that will be discussed below, a 1 is a smoothness 
indicator that is close to zero in smooth regions and close to 0.5 near disconti- 
nuities. Note that we distinguish the size of the total time step: 

At = t n+1 -t n , (95) 

from the time step of a particular substep of the Strang splitting: 

At = t™ +1 - r n . (96) 



22 



We add this additional diffusion term to each weakly hyperbolic equation in our 
dimensionally split algorithm in the following way: 



Sub-problem 1: 



A 4. U i\ ^, 11 J\ ^, — £ A ~ 



A 3 



u x A 2 x = 0, 
wM 3 , = 0, 



(97) 



Sub-problem 2: A\ + u 2 'A) y = 0, 



J\ a. U J\ U J± — £ J± 

; L iy it/ ;U 

A 3 t + u 2 A% = 0, 



(98) 



Sub-problem 3: 



4 

4 
A\ 



l A\ = 0, 



u 3 A% 
u x A\ 



= 0, 

u 2 A\ 



(99) 



3 A3 



£ 6 A 



Note that we diffuse only in the direction of the dimensionally split solve. 

Consider again only the discretization used in Sub-problem (97). The nu- 
merical update in the weakly hyperbolic solve can be written as 



where [.A 1 ] ... is given from update (91). Stability requires that 

At 1 

0< 2 — va< -, 
~At ~ 2 

which can be achieved if we assume that At < At and we take 



< a < 



1 



and 



0< v< 



1 



(100) 



(101) 



(102) 



In our proposed method, a is a smoothness indicator and v is a user-defined pa- 
rameter (typically i/<l) that can be tuned to control the amount of numerical 
dissipation across discontinuities in the magnetic field. 

The smoothness indicator a is computed according to the formulas 



max 



«i 



Ol 



a\ + a r 



with 



«1 



— <,e+[A i: j k ^4i_ijfc) > , a T — ie + [A i+1 j k A i:jk ) > 



(103) 



(104) 



The parameter e is introduced only to avoid division theory and is in practice 
taken to be e = 10~ 8 . The smoothness indicator is designed to distinguish 
between the following cases: 
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1. If A 1 is smooth for x € (xi_i,Xj+i), then 

A \ik - A l-ijk = 0(Ax) 7 and A] +ljk - A] jk = O(Ax). 
In this case one can show that 

.1' 



lim lim a 

Az->0 e-5-0 
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Ax = O(Ax), 



A \ 

which shows that the overall numerical method still retains 0(At 2 , Ax 2 ) 
accuracy. 

2. If A 1 is non-smooth in (x,_i,Xj) and smooth in (x,,x,+i), then a r 3> a\ 
and a w g • 

3. If A 1 is non-smooth in (xj,X,+i) and smooth in (xj_i,x,), then ai ^> a r 
and a « \ . 

In all cases a < 1/2, which guarantees that the numerical update will be stable 
up to CFL number one. 

5.4- Discretization of the 2. 5-dimensional problem 

In order to construct and analyze methods for the three-dimensional MHD 
equations, it is instructive to also look at the so-called 2. 5-dimensional case. 
In this case we assume that u and B are three-dimensional vectors, but all 
conserved quantities [p, pu, £, B] T are functions of only two spatial variables, 
say x= {x,y) T . 

Since B\ — 0, we can obtain a divergence free magnetic field by employing a 
constrained transport algorithm that only updates B 1 and B 2 , and which treats 
B 3 as just another conserved variable in the base scheme. In other words, one 
could simply update the transport equation for A 3 , 

A 3 t + u l A 3 x + u 2 A 3 y = 0, (105) 

via the method of [35] , then compute the first two components of the magnetic 
field from discrete analogs of B 1 — A 3 , and B 2 = — A 3 , and let the base scheme 
update B since it doesn't enter into the 2.5D divergence constraint: 

B% + B 2 y = 0. 

An alternative approach to handling the 2.5D case, one that will allow us to 
test the important features of the proposed 3D scheme, is to consider the full 
magnetic vector potential evolution equation in the xy-planc: 

(106) 
(107) 
(108) 



) t - u 2 A 2 x - u 3 A 3 x + u 2 A] y 


= 


2 t +u 1 A 2 x -u 1 A] y -u 3 A 3 y 


= 


i4f t + «Uf x + « 2 Af tf 


= 
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Approximate solutions to these equations can be computed via the method 
outlined in §5.2 and §5.3. The magnetic field can be obtained via discrete 
analogs of 

B 1 ^A 3 y , B 2 = -A%, and B 3 = A 2 X - A) y . 

Because the numerical update based on the scheme presented in §5.2 and §5.3 is 
not equivalent to the method of [35] in 2.5D, we expect to see some differences 
in the solution of the magnetic field. However, since the two methods are very 
closely related, we expect to obtain similar results. Numerical examples for 
the 2.5D case comparing the 2D constrained transport method of [35] and the 
proposed 3D method are discussed in §6.1.1 and §6.1.2. 

6. Numerical experiments 

We present several numerical experiments in this section. All the numerical 
tests are carried out in the CLAWPACK software package [22] , which can be freely 
downloaded from the web. In particular, the current work has been incorporated 
into the MHDCLAW extension of CLAWPACK, which was originally developed by 
Rossmanith [34]. This extension can also be freely downloaded from the web. 

6.1. Test cases in 2.5D 

We begin by presenting numerical results for two problems in 2.5 dimensions 
(i.e., the solution depends only on the independent variables (x,y,t), but all 
three components of the velocity and magnetic field vectors are non-zero) . The 
two examples that are considered are: 

1. Smooth Alfven wave; and 

2. Cloud-shock interaction. 

The first example involves an infinitely smooth exact solution of the ideal MHD 
system. The second example involves the interaction of a high-pressure shock 
with a high-density bubble (i.e., a discontinuous example). 

6.1.1. Smooth Alfven wave problem, 

We first verify the order of convergence of the constrained transport method 
for a smooth circular polarized Alfven wave that propagates in direction n = 
(cos (f>, sin <j>, 0) T towards the origin. This problem has been considered by several 
authors (e.g., [35, 40]), and is a special case (9 = 0) of the 3D problem described 
in detail in §6.2.1. Here we take <j> = tan _1 (0.5) and solve on the domain: 

(x,y)€ 0,^-r x 0,^- . (109) 

COS0J [ Sllllp 

The solution consists of a sinusoidal wave propagating at constant speed without 
changing shape, thus making it a prime candidate to verify order of accuracy. 

In Table 1, we show the convergence rates for the 2.5D test problem. Here 
all three components of the magnetic field have been updated in the constrained 



2o 



Mesh 


Loo Error in B 1 


Loo Error in B 2 


Loo Error in B 3 


64 x 128 


6.778 x 1(T 4 


2.393 x 10~ 3 


1.284 x 10~ 2 


128 x 256 


1.690 x 10~ 4 


5.969 x 10~ 4 


3.203 x 10~ 3 


256 x 512 


4.221 x 10~ 5 


1.492 x 10~ 4 


8.004 x 10~ 4 


512 x 1024 


1.057 x 10~ 5 


3.741 x 10~ 5 


2.000 x 10~ 4 


Order 


2.001 


2.003 


2.005 



Mesh 


Loo Error in A 1 


Loo Error in A 2 


Loo Error in A 3 


64 x 128 


1.302 x 10~ 2 


1.288 x 10^2 


1.453 x 10~ 2 


128 x 256 


3.260 x 10~ 3 


3.217 x 10~ 3 


3.633 x 10~ 3 


256 x 512 


8.213 x 10~ 4 


8.025 x 10~ 4 


9.081 x 10~ 4 


512 x 1024 


2.089 x 10~ 4 


1.997 x 10~ 4 


2.270 x 10~ 4 


Order 


1.987 


2.004 


2.001 



Table 1: Error tables for the 2.5D Alfven problem at time t = 1.5 using proposed dimensional 
split scheme. This table shows that all components of the magnetic field and all components 
of the magnetic potential converge at second order accuracy. 



transport method and the magnetic potential was approximated using a dimen- 
sional split method for (106)-(108). The table clearly shows that the proposed 
method is second-order accurate in all of the magnetic field components, as well 
as the magnetic potential components. The reported errors are of the same 
magnitude as those reported in Table 7.1 of [35] (page 1786). 

6.1.2. Cloud-shock interaction problem, 

Next we consider a problem with a strong shock interacting with a relatively 
large density jump in the form of a high-density bubble that is at rest with its 
background before the shock-interaction. This problem has also been considered 
by several authors (e.g., [35, 40]), however, in previous work a solution was 
always obtained by treating B 3 as a standard conserved variable. Here we 
compare such an approach in the form of the method described in [35] , against 
the method that is proposed in the current work. In our new method all three 
components of the magnetic field are computed from derivatives of the magnetic 
vector potential. The details of the initial conditions are written in §6.3.2. 

In Figure 1 we show numerical results for the 2.5 dimensional problem. 
For the newly proposed method we take v = 0.02 as our artificial diffusion 
parameter. We compare the scheme from [35] that only updates B 1 and B 2 
using the scalar evolution equation for the vector potential (105) with the new 
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Figure 1: The 2.5D cloud-shock interaction problem. Shown here are the out-of-plane mag- 
netic field at time t = 0.06 as calculated on a 512x512 mesh using (a) the scheme of [35] that 
only uses A 3 , and (b) the proposed scheme using the full vector potential A. 



constrained transport method that updates all components of the magnetic field. 
Although these methods compute B 3 in very different ways, using very different 
limiters, these plots show that they nonetheless produce similar solutions. This 
result gives us confidence that the proposed method is able to accurately resolve 
strong shocks despite the somewhat unusual approach for approximating the 
magnetic field. 

6.2. Test cases in 3D 

We present numerical results for three-dimensional versions of three classical 
MHD test problems: 

1. Smooth Alfven wave; 

2. Rotated shock tube problem; 

3. Orszag-Tang vortex; and 

4. Cloud-shock interaction. 

Two-dimensional versions of these problems have been studied by many authors, 
see for instance [35, 40]. 

Our test calculations will all be based on the splitting (77)-(79). We also 
carried out several tests with (74)-(76); however, we will not report those here. 
The two methods produced comparable results. 

6.2.1. Smooth Alfven wave problem 

We first verify the order of convergence of the constrained transport method 
for a smooth circular polarized Alfven wave that propagates in direction n = 
(cos c/> cos 9 , sin <p cos 9 , sm 9) T towards the origin. Here <fi is an angle with re- 
spect to the x-axis in the a;y-plane and 9 is an angle with respect to the xz- 
plane. Initial values for the velocity and the magnetic field are specified in the 
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direction n as well as the orthonormal directions t = (— sin</>, cos<j>, 0) T and 
r = (— cos <j) sin 9, — sin <\> sin 9, cos 9) T : 



where 



u(0, x) = u n n + u l t + u r r, 
B(0,x) = B"n + S t t + S r r, 



u* =B* = 0.1sin(27rn-x), 
M r = B r = 0.1cos(27rn-x). 



The initial density and pressure are constant and set to 
p(0,x) = l and p(0,x)=0.1, 



(110) 
(111) 



(112) 
(113) 
(114) 



(115) 



respectively. This choice guarantees that the Alfven wave speed is \va\ = 
Bn/y/p — 1) which means that the flow agrees with the initial state whenever 
the time is an integer value. The computational domain is taken to be 



Q 



0. 



('( IS COS f 



0. 



0,-7 



sin^ 



(116) 



with periodic boundary conditions imposed on the conserved variables in all 
three coordinate directions: (p, pu, £ , B). 

Our initial condition for the magnetic potential is 



A 1 ^,*.) = z(B 2 ) sin0sin(27rn-x), 

A 2 (Q,x) = x(B 3 ) H cos0sin(27rn-x), 



A 3 (0,x) = y(B 1 ) + 



1 



■ cos(2-7rn • x), 



(117) 
(118) 
(119) 



20tt cos 9 
diere (B) denotes the average magnetic field over the computational domain 



1 

W\ 



B(i, x, y, z) dx dydz — (cos <f> cos 9, sin <f) cos 9, sin 9) . (120) 



We note that even though the magnetic field is time-dependent, its average, 
(B), is time- independent. Therefore, the magnetic potential consists of a linear 
(time-independent) and a periodic (time-dependent) part. Boundary conditions 
for the magnetic potential are handled by applying periodic boundary condi- 
tions on the periodic part and linear extrapolation for the linear part (linear 
extrapolation is exact in this case). 

In Table 2 we show the results of a numerical convergence study in the three- 
dimensional case with 6 = <fi = tan -1 (0.5) w 26.5651°. These results confirm 
that our method is second order accurate. 
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Mesh 


Loo Error in B 1 


Loo Error in B 2 


Loo Error in B 3 


16 x 32 x 32 


1.022 x 1CT 2 


2.787 x 10~ 2 


2.382 x 10~ 2 


32 x 64 x 64 


2.577 x 1(T 3 


7.075 x 10~ 3 


6.101 x 10~ 3 


64 x 128 x 128 


6.487 x 1(T 4 


1.782 x 10~ 3 


1.549 x 10~ 3 


Order 


1.990 


1.989 


1.978 



Mesh 


Loo Error in A 1 


Loo Error in A 2 


Loo Error in A 3 


16 x 32 x 32 


1.902 x lO" 1 


1.460 x 10" 1 


1.187 x lO" 1 


32 x 64 x 64 


4.816 x 10~ 2 


3.747 x 10- 2 


2.995 x 10~ 2 


64 x 128 x 128 


1.220 x 10~ 2 


9.528 x 10~ 3 


7.564 x 10~ 3 


Order 


1.981 


1.975 


1.985 



Tabic 2: Error tables for the 3D Alfvcn problem at time t = 1. This table shows that all 
components of the magnetic field and all components of the magnetic potential converge at 
second order accuracy. The approximate order of accuracy is computed by taking the log 2 of 
the ratio of the error on the 32 X 64 X 64 mesh to the error on the 64 X 128 X 128 mesh. 



6.3. Rotated shock tube problem 

We now consider the shock tube problem described in [27, 38]. In this test 
case a one-dimensional shock tube problem is computed in a three-dimensional 
domain. In order to define the problem, we consider a coordinate transformation 
of the form 



cos(a) cos(/3) cos(a) sin(/3) sin(a) 

- sin(/3) cos(/3) 

— sin(a) cos(/3) — sin(a) sin(/3) cos(a) 



(121) 



which maps the Cartesian coordinates (x, y, z) T to the rotated coordinate system 
(£,7VX) T - Vectors transform as follows 



B x 

B» 
B z 



cos(a) cos(/3) — sin(/3) — sin(a) cos(/3) 

cos(a) sin(/3) cos(/3) — sin(a) sin(/3) 

sin(a) cos(a) 



B'i 
B< 



(122) 



The Ricmann initial data is taken to be 



(P. 



M s , U ' U 



,p,Be,B*,B<)(0,x) 



1.08, 1.2, 0.01, 0.5, 0.95 
> 1, o, 0, 0, 1 



2 

/In' 



2 3.6 

"StT ' \/iTT ' 

4 2 

Iff ' \/4¥ . 



2 

/4tt 



if e<o, 
if e>o, 



(123) 
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where the velocity vector and the magnetic field are given in the rotated coor- 
dinate frame. As the initial condition for the magnetic potential (in the rotated 
coordinate frame) we use 

(A«, A\ A<) (0, x) = (0, £B<, 7,13* - ££") . (124) 

The solution of the Riemann problem at later times remains a function 
of £. This fact is used to extend the computed values to ghost cell values. 
The computational domain is Q = [-0.75,0.75] x [0,0.015625] x [0,0.015625], 
which is discretized using 768 x 8 x 8 grid cells. By using a = tan -1 (0.5) 
and j3 = tan -1 (0.25 cos a), we obtain simple formulas for the extension of the 
numerical solution to the ghost cells. To define ghost cell values in the Cartesian 
y-direction, we can for instance use the relation q(i,j, k) = q(i + l,j — 2, k). In 
the z-direction we can use q(i,j, k) = q(i + l,j, k — 4). The ghost cell values in 
the x-direction were computed via extrapolation. We present numerical results 
at time t = 0.2 as scatter plots of the three-dimensional solution plotted as a 
function of £. 

As already documented by other authors, we also observe some oscillations 
in particular in the B* component of the magnetic field. These oscillations do 
not appear if the shock tube initial data are aligned with the mesh. In Figure 
2 we show approximations of the solution at time t = 0.2. These plots are 
scatterplots, where the output of the three-dimensional computation is plotted 
as a function of £. The solid line in these plots is obtained by computing 
solutions of the one-dimensional Riemann problem on a very fine mesh (using 
10 grid cells). The computation was performed using the minmod limiter and 
the diffusive limiter with v = 0.05. As expected, with a less diffusive limiter the 
amplitude of the oscillations becomes stronger, by increasing v, the amplitude 
of the oscillations becomes weaker. 

6.3.1. Orszag-Tang vortex 

Our next test problem is the Orszag-Tang vortex problem. This is a stan- 
dard test case for the two-dimensional MHD equations. Here we add a small 
perturbation to the initial velocity field that depends on the vertical direction. 

The initial condition is 

p(0,x)= 7 2 , p(0,x)= 7) (125) 

u(0,x) = ( — (1 + esinz) siny, (1 + esinz) sinx, esin(z)) , (126) 



B(0,x)= (-amy, sin(2x), o) . (127) 

Here we used e = 0.2. The initial condition for the magnetic potential is 

A(0,x) = (O, 0, cosy + cos(2a;)) . (128) 

The computational domain is a cube with side length 2ir. Periodicity is imposed 
in all three directions. 
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In Figures 3-6 we show a sequence of schlieren plots of the pressure for several 
slices of data in the rry-planc for z = n/2, ir, 3ir/2. Those computations have 
been performed with the high-resolution constrained transport method with 
monotonized central limiter and the diffusive limiter using v — 0.05. In Figures 
7, 8 we show schlieren plots of pressure in the xy-plane for z — n at different 
times. 

6.3.2. Cloud-shock interaction problem 

Finally we consider the cloud-shock interaction problem. The initial condi- 
tions consist of a shock that is initially located at x — 0.05: 

(p,u 1 ,u 2 ,u 3 ,p,B 1 ,B 2 ,B 3 )(0,x) 

(3.86859, 11.2536, 0, 0, 167.345, 0, 2.1826182, -2.1826182) if x < 0.05, 
(1, 0, 0, 0, 1, 0, 0.56418958, 0.56418958) if i> 0.05, 

(129) 

and a spherical cloud of density p = 10 with radius r = 0.15 and centered at 
(0.25,0.5,0.5). The cloud is in hydrostatic equilibrium with the fluid to the 
right of the shock. The initial condition for the magnetic potential is 

, . _ f (2.1826182 y, 0, -2.1826182 (a; - 0.05)) T if i< 0.05, 
^ ' X ^ ~ | (-0.56418958 y, 0, 0.56418956 (x - 0.05)) T if i> 0.05. 

(130) 

This test problem can be computed on the unit cube with inflow boundary 
conditions on the left side and outflow conditions on all other sides. Instead we 
make use of the symmetry of the problem and compute the problem only in a 
quarter of the full domain, i.e. in [0, 1] x [0.5, 1] x [0.5, 1] and impose reflecting 
boundary conditions at the lower boundary in the y and z-dircctions. For the 
conserved quantities, the reflecting boundary condition is implemented in the 
standard way, i.e., by copying the values of the conserved quantities from the 
first grid cells of the flow domain to the ghost cells. The normal momentum 
component in the ghost cells is negated. The components of the magnetic 
potential are linearly extrapolated to the ghost cells. 

In Figures 9 and 10, we show a sequence of schlieren plots of the density in 
the xy-plane for z — 0.5 and the xz-plane for y = 0.5. The three-dimensional 
radial symmetric solution structure compares well with previously shown two- 
dimensional computations. In Figures 11 and 12 we show the schlieren plots of 
density in the xz-plane. Here the diffusive limiter described in §5.2 was used 
with v = 0.02. 



7. Conclusions 

We have presented a new constrained transport method for the three-dimen- 
sional MHD equations. Depending on the gauge condition, we discussed differ- 
ent possible evolution equations for the magnetic potential. All of these gauge 
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conditions implicate some difficulties for the discretization. The condition used 
here leads to a weakly hyperbolic system for the evolution of the magnetic po- 
tential. We discretized this system with a splitting approach. For the MHD 
equations we used a wave propagation method. Several numerical tests confirm 
the robustness and accuracy of the resulting constrained transport scheme. 

Our method is fully explicit, as well as fully unstaggered, and therefore well- 
suited for adaptive mesh refinement and parallelization. These generalizations 
will be the focus of future work. 
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Figure 2: The rotated Riemann problem. Scatterplots of various components of solution at 
time t = 0.2 using the constrained transport algorithm. The solid line in each panel is a 
highly-resolved ID simulation. 
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Pressure at time t=1 




Figure 3: The 3D Orszag-Tang problem. Shown in this figure are schlieren slices of pressure 
at time t = 1 and at various z values (z = n/2, n, 37r/2). This solution was computed on a 
mesh with 150 X 150 X 150 grid cells. 
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Pressure at time t=2 




Figure 4: The 3D Orszag-Tang problem. Shown in this figure are schlieren slices of pressure 
at time t = 2 and at various z values (z = n/2, n, 37r/2). This solution was computed on a 
mesh with 150 X 150 X 150 grid cells. 
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Pressure at time t=3 
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Figure 5: The 3D Orszag-Tang problem. Shown in this figure are schlieren slices of pressure 
at time t = 3 and at various z values (z = n/2, n, 37r/2). This solution was computed on a 
mesh with 150 X 150 X 150 grid cells. 
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Pressure attimet=3.5 




Figure 6: The 3D Orszag-Tang problem. Shown in this figure are schlieren slices of pressure 
at time t = 3.5 and at various z values (z = 7r/2, n, 37r/2). This solution was computed on a 
mesh with 150 X 150 X 150 grid cells. 
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Pressure at time t=1 
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Figure 7: The 3D Orszag-Tang problem. Schliercn plots of pressure at time (a) t = 1 and (b) 
t = 2 computed on a mesh with 150 X 150 X 150 grid cells. Slices of the solution at z = n in 
the xy-planc arc shown. 
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Pressure at time t=3 
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Figure 8: The 3D Orszag-Tang vortex problem. Shown in these panels are schlieren plots of 
pressure in the xj/-planc at (a) time t = 3 and z = it and (b) t = 3.5 and z = n. These 
solutions were computed on a mesh with 150 X 150 X 150 grid cells. 
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Figure 9: The 3D cloud-shock interaction problem. Schlicrcn plots of density for the problem 
using a mesh with 200 X 100 X 100 grid cells at time (a) t = and (b) t = 0.02. Shown here 
is the solution in two selected orthogonal planes. 
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Figure 10: The 3D cloud-shock interaction problem. Schlieren plots of density using a mesh 
with 200 X 100 X 100 grid cells at time (a) t = 0.04 and (b) t = 0.06. Shown here is the 
solution in two selected orthogonal planes. 
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Figure 11: The 3D cloud-shock interaction problem. Sequence of schlieren plots of density in 
the 22-plane for y = 0.5 at time (a) t = and (b) t = 0.02. 
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Density at time t=0.04 
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Figure 12: The 3D cloud-shock interaction problem. Sequence of schlieren plots of density in 
the iz-plane for y = 0.5 at time (a) t = 0.04 and (b) t = 0.06. 
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